##################
# Plot Estimates #
##################  
pdf(paste("Estimates",str.sample,str.rho,str.cb,str.robust,str.boot,".pdf",sep=""));

par.current <- par();

# 1. Counterfactual level 1 #

# Main Parameters
par(mfrow=c(3,2),mar=c(5,4,4,2),oma=c(0,0,2,0));  
for(i in 1:length(para.namelist)){
  para.name <- para.namelist[i];
  matplot(q.plot/100,
          cbind(-lbound.counter1[qn.plot,para.name],
                -theta.counter1[qn.plot,para.name],
                -ubound.counter1[qn.plot,para.name]),
          col=c("red","steelblue4","red"), type="l", lty=1, ylim=ylimlist[i,],
          xlab="Quantile Index", ylab="", main=para.name);
  abline(h=0);
  abline(h=beta.Heckman.counter1[para.name]/sigma.Heckman.counter1,
         lty=2);
  if(i%%6==0 | i==length(para.namelist)){
    #mtext(paste("Estimates of Parameters, ", main.str.counter1, sep=""), outer=T, cex=1.2);
    mtext(main.str.counter1, outer=T, cex=1.2);};
};



# Rho #
if(rho.spec==1){
  par(mfrow=c(1,1),mar=c(5,4,4,2),oma=c(0,0,2,0));
  matplot(q.plot/100,
          cbind(-tanh(theta.counter1[qn.plot,"delta.(Intercept)"])-t(crt.counter1*sqrt(var.counter1))[qn.plot,"delta.(Intercept)"]/(cosh(theta.counter1[qn.plot,"delta.(Intercept)"]))^2,
                -tanh(theta.counter1[qn.plot,"delta.(Intercept)"]),
                -tanh(theta.counter1[qn.plot,"delta.(Intercept)"])+t(crt.counter1*sqrt(var.counter1))[qn.plot,"delta.(Intercept)"]/(cosh(theta.counter1[qn.plot,"delta.(Intercept)"]))^2),
          col=c("red","steelblue4","red"), type="l", lty=1, ylim=c(-1,1), 
          xlab="Quantile Index", ylab="");
  abline(h=0);
  abline(h=rho.Heckman.counter1,lty=2);   
  #mtext(paste("Estimates of Unobserved Selection, ", main.str.counter1, sep=""), outer=T, cex=1.2);
  mtext(main.str.counter1, outer=T, cex=1.2);
  
}else if(rho.spec==2 | rho.spec==3 | rho.spec==9 | rho.spec==10 | rho.spec==12 | rho.spec==14 | rho.spec==19){
  par(mfrow=mfrow.current,mar=c(5,4,4,2),oma=c(0,0,2,0)); 
  for(i in 1:length(delta.each)){
    delta.para <- delta.each[i];
    delta.name <- paste("delta.",delta.para,sep="");
    matplot(q.plot/100,
            cbind(-tanh(theta.counter1[qn.plot,delta.name])-t(crt.counter1*sqrt(var.counter1))[qn.plot,delta.name]/(cosh(theta.counter1[qn.plot,delta.name]))^2,
                  -tanh(theta.counter1[qn.plot,delta.name]),
                  -tanh(theta.counter1[qn.plot,delta.name])+t(crt.counter1*sqrt(var.counter1))[qn.plot,delta.name]/(cosh(theta.counter1[qn.plot,delta.name]))^2),
            col=c("red","steelblue4","red"), type="l", lty=1, ylim=c(-1,1),
            xlab="Quantile Index", ylab="", main=bquote(rho(y) ~ "for" ~ .(delta.para)));
    abline(h=0);
    abline(h=rho.Heckman.counter1,lty=2);   
    if(rho.spec==2 & i%%6==0){#mtext(paste("Estimates of Unobserved Selection, ", main.str.counter1, sep=""), outer=T, cex=1.2)
      mtext(main.str.counter1, outer=T, cex=1.2)};
  };
  if(rho.spec!=2){
    #mtext(paste("Estimates of Unobserved Selection, ", main.str.counter1, sep=""), outer=T, cex=1.2)
    mtext(main.str.counter1, outer=T, cex=1.2)};
}else if(rho.spec==5 | rho.spec==6 | rho.spec==7 | rho.spec==8 | rho.spec==11 | rho.spec==13 | rho.spec==17 | rho.spec==20 | rho.spec==21){
  par(mfrow=mfrow.current,mar=c(5,4,4,2),oma=c(0,0,2,0));
  for(i in 1:length(delta.namelist)){
    delta.para <- delta.namelist[i];   
    matplot(q.plot/100,
            cbind(-lbound.counter1[qn.plot,paste("delta.",delta.para,sep="")],
                  -theta.counter1[qn.plot,paste("delta.",delta.para,sep="")],
                  -ubound.counter1[qn.plot,paste("delta.",delta.para,sep="")]),
            col=c("red","steelblue4","red"), type="l", lty=1, ylim=ylim.current,
            xlab="Quantile Index", ylab="", 
            main=bquote("Effect"~"of"~.(delta.para)~"on"~delta(y)));
    abline(h=0);
  };   
  #mtext(paste("Estimates of Unobserved Selection, ", main.str.counter1, sep=""), outer=T, cex=1.2);
  mtext(main.str.counter1, outer=T, cex=1.2);
}else if(rho.spec==4){
  par(mfrow=mfrow.current,mar=c(5,4,4,2),oma=c(0,0,2,0)); 
  for(i in 1:length(delta.each)){
    delta.para <- delta.each[i];
    delta.name <- paste("delta.",delta.para,sep="");
    matplot(q.plot/100,
            cbind(-tanh(theta.counter1[qn.plot,delta.name])-t(crt.counter1*sqrt(var.counter1))[qn.plot,delta.name]/(cosh(theta.counter1[qn.plot,delta.name]))^2,
                  -tanh(theta.counter1[qn.plot,delta.name]),
                  -tanh(theta.counter1[qn.plot,delta.name])+t(crt.counter1*sqrt(var.counter1))[qn.plot,delta.name]/(cosh(theta.counter1[qn.plot,delta.name]))^2),
            col=c("red","steelblue4","red"), type="l", lty=1, ylim=c(-1,1),
            xlab="Quantile Index", ylab="", main=bquote(rho(y) ~ "for" ~ .(delta.para)));
    abline(h=0);
    abline(h=rho.Heckman.counter1,lty=2);   
    if(i%%6==0){#mtext(paste("Estimates of Unobserved Selection, ", main.str.counter1, sep=""), outer=T, cex=1.2)
      mtext(main.str.counter1, outer=T, cex=1.2)};
  };
  if(length(delta.each)%%6!=0){#mtext(paste("Estimates of Unobserved Selection, ", main.str.counter1, sep=""), outer=T, cex=1.2)
    mtext(main.str.counter1, outer=T, cex=1.2)};
  
  par(mfrow=c(1,1),mar=c(5,4,4,2),oma=c(0,0,0,0));
  mycols <- topo.colors(length(rho.name));
  matplot(q.plot/100, -tanh(theta.counter1[qn.plot, rho.name]),
          col=mycols, type="l", lty=1, ylim=c(-1,1),
          xlab="Quantile Index", ylab="", main=main.str.counter1);
  #paste("Estimates of Unobserved Selection by Years, ", main.str.counter1, sep="")
  abline(h=0);
  image.plot(legend.only=T, zlim=c(1978,2013), col=mycols, smallplot=c(0.8,0.82,0.2,0.4),
             axis.args=list(at=c(1978,1995,2013),labels=c("1978","1995","2013")));
}else if(rho.spec==15){
  par(mfrow=c(1,1),mar=c(5,4,4,2),oma=c(0,0,0,0));
  allt <- unique(mysample$year_res);
  mycols <- topo.colors(length(allt));
  matplot(q.plot/100,
          -tanh(replicate(length(allt), theta.counter1[qn.plot,"delta.(Intercept)"]) + 
                  theta.counter1[qn.plot,"delta.year_res"] %o% allt),
          col=mycols, type="l", lty=1, ylim=c(-1,1),
          xlab="Quantile Index", ylab="", main=main.str.counter1);
  #paste("Estimates of Unobserved Selection by Years, ", main.str.counter1, sep="")
  abline(h=0);
  image.plot(legend.only=T, zlim=c(1978,2013), col=mycols, smallplot=c(0.8,0.82,0.2,0.4),
             axis.args=list(at=c(1978,1995,2013),labels=c("1978","1995","2013")));
  
  par(mfrow=mfrow.current,mar=c(5,4,4,2),oma=c(0,0,2,0));
  for(i in 1:length(delta.namelist)){
    delta.para <- delta.namelist[i];   
    matplot(q.plot/100,
            cbind(-lbound.counter1[qn.plot,paste("delta.",delta.para,sep="")],
                  -theta.counter1[qn.plot,paste("delta.",delta.para,sep="")],
                  -ubound.counter1[qn.plot,paste("delta.",delta.para,sep="")]),
            col=c("red","steelblue4","red"), type="l", lty=1, 
            xlab="Quantile Index", ylab="", 
            main=bquote("Effect"~"of"~.(delta.para)~"on"~delta(y)));
    abline(h=0);
  }; 
  #mtext(paste("Estimates of Unobserved Selection, ", main.str.counter1, sep=""), outer=T, cex=1.2);
  mtext(main.str.counter1, outer=T, cex=1.2);
}else if(rho.spec==16){
  par(mfrow=c(1,1),mar=c(5,4,4,2),oma=c(0,0,0,0)); 
  allt <- unique(mysample$year_res);
  mycols <- topo.colors(length(allt));
  matplot(q.plot/100,
          -tanh(replicate(length(allt), theta.counter1[qn.plot,"delta.(Intercept)"]) + 
                  theta.counter1[qn.plot,"delta.year_res"] %o% allt +
                  theta.counter1[qn.plot,"delta.I(year_res^2)"] %o% allt^2),
          col=mycols, type="l", lty=1, ylim=c(-1,1),
          xlab="Quantile Index", ylab="", main=main.str.counter1);
  #paste("Estimates of Unobserved Selection by Years, ", main.str.counter1, sep="")
  abline(h=0);
  image.plot(legend.only=T, zlim=c(1978,2013), col=mycols, smallplot=c(0.8,0.82,0.2,0.4),
             axis.args=list(at=c(1978,1995,2013),labels=c("1978","1995","2013")));
  
  par(mfrow=mfrow.current,mar=c(5,4,4,2),oma=c(0,0,2,0));
  for(i in 1:length(delta.namelist)){
    delta.para <- delta.namelist[i];   
    matplot(q.plot/100,
            cbind(-lbound.counter1[qn.plot,paste("delta.",delta.para,sep="")],
                  -theta.counter1[qn.plot,paste("delta.",delta.para,sep="")],
                  -ubound.counter1[qn.plot,paste("delta.",delta.para,sep="")]),
            col=c("red","steelblue4","red"), type="l", lty=1, 
            xlab="Quantile Index", ylab="", 
            main=bquote("Effect"~"of"~.(delta.para)~"on"~delta(y)));
    abline(h=0);
  }; 
  #mtext(paste("Estimates of Unobserved Selection, ", main.str.counter1, sep=""), outer=T, cex=1.2);
  mtext(main.str.counter1, outer=T, cex=1.2);
}else if(rho.spec==18){
  par(mfrow=c(2,1),mar=c(5,4,4,2),oma=c(0,0,2,0)); 
  allt <- unique(mysample$year_res);
  mycols <- topo.colors(length(allt));
  matplot(q.plot/100,
          -tanh(replicate(length(allt), theta.counter1[qn.plot,"delta.(Intercept)"]) + 
                  theta.counter1[qn.plot,"delta.year_res"] %o% allt),
          col=mycols, type="l", lty=1, ylim=c(-1,1),
          xlab="Quantile Index", ylab="", main=bquote(rho(y)~"for Single"));
  abline(h=0);
  #mtext(paste("Estimates of Unobserved Selection by Years, ", main.str.counter1, sep=""), outer=T, cex=1.2);
  mtext(main.str.counter1, outer=T, cex=1.2);
  
  matplot(q.plot/100,
          -tanh(replicate(length(allt), theta.counter1[qn.plot,"delta.(Intercept)"] + theta.counter1[qn.plot,"delta.couple"]) + 
                  (theta.counter1[qn.plot,"delta.year_res"] + theta.counter1[qn.plot,"delta.year_res:couple"]) %o% allt),
          col=mycols, type="l", lty=1, ylim=c(-1,1),
          xlab="Quantile Index", ylab="", main=bquote(rho(y)~"for Married"));
  abline(h=0);
  image.plot(legend.only=T, zlim=c(1978,2013), col=mycols, smallplot=c(0.8,0.82,0.35,0.45),
             axis.args=list(at=c(1978,2013),labels=c("1978","2013")));  
  
  par(mfrow=mfrow.current,mar=c(5,4,4,2),oma=c(0,0,2,0));
  for(i in 1:length(delta.namelist)){
    delta.para <- delta.namelist[i];   
    matplot(q.plot/100,
            cbind(-lbound.counter1[qn.plot,paste("delta.",delta.para,sep="")],
                  -theta.counter1[qn.plot,paste("delta.",delta.para,sep="")],
                  -ubound.counter1[qn.plot,paste("delta.",delta.para,sep="")]),
            col=c("red","steelblue4","red"), type="l", lty=1, 
            xlab="Quantile Index", ylab="", 
            main=bquote("Effect"~"of"~.(delta.para)~"on"~delta(y)));
    abline(h=0);
  }; 
  #mtext(paste("Estimates of Unobserved Selection, ", main.str.counter1, sep=""), outer=T, cex=1.2);
  mtext(main.str.counter1, outer=T, cex=1.2);
};


# 2. Counterfactual level 2 #

# Main Parameters #
par(mfrow=c(3,2),mar=c(5,4,4,2),oma=c(0,0,2,0)); 
for(i in 1:length(para.namelist)){
  para.name <- para.namelist[i];
  matplot(q.plot/100,
          cbind(-lbound.counter2[qn.plot,para.name],
                -theta.counter2[qn.plot,para.name],
                -ubound.counter2[qn.plot,para.name]),
          col=c("red","steelblue4","red"), type="l", lty=1, ylim=ylimlist[i,],
          xlab="Quantile Index",ylab="", main=para.name);
  abline(h=0);
  abline(h=beta.Heckman.counter2[para.name]/sigma.Heckman.counter2, lty=2);
  if(i%%6==0 | i==length(para.namelist)){#mtext(paste("Estimates of Parameters, ", main.str.counter2, sep=""), outer=T, cex=1.2);
    mtext(main.str.counter2, outer=T, cex=1.2);};
};



# Rho #
if(rho.spec==1){
  par(mfrow=c(1,1),mar=c(5,4,4,2),oma=c(0,0,2,0));
  matplot(q.plot/100,
          cbind(-tanh(theta.counter2[qn.plot,"delta.(Intercept)"])+t(crt.counter2*sqrt(var.counter2))[qn.plot,"delta.(Intercept)"]/(cosh(theta.counter2[qn.plot,"delta.(Intercept)"]))^2,
                -tanh(theta.counter2[qn.plot,"delta.(Intercept)"]),
                -tanh(theta.counter2[qn.plot,"delta.(Intercept)"])-t(crt.counter2*sqrt(var.counter2))[qn.plot,"delta.(Intercept)"]/(cosh(theta.counter2[qn.plot,"delta.(Intercept)"]))^2),
          col=c("red","steelblue4","red"), type="l", lty=1, ylim=c(-1,1), xlab="Quantile Index",ylab="");
  abline(h=0);
  abline(h=rho.Heckman.counter2, lty=2);
  #mtext(paste("Estimates of Unobserved Selection, ", main.str.counter2, sep=""), outer=T, cex=1.2);
  mtext(main.str.counter2, outer=T, cex=1.2);
  
}else if(rho.spec==2 | rho.spec==3 | rho.spec==9 | rho.spec==10 | rho.spec==12 | rho.spec==14 | rho.spec==19){
  par(mfrow=mfrow.current,mar=c(5,4,4,2),oma=c(0,0,2,0)); 
  for(i in 1:length(delta.each)){
    delta.para <- delta.each[i];
    delta.name <- paste("delta.",delta.para,sep="");
    matplot(q.plot/100,
            cbind(-tanh(theta.counter2[qn.plot,delta.name])-t(crt.counter2*sqrt(var.counter2))[qn.plot,delta.name]/(cosh(theta.counter2[qn.plot,delta.name]))^2,
                  -tanh(theta.counter2[qn.plot,delta.name]),
                  -tanh(theta.counter2[qn.plot,delta.name])+t(crt.counter2*sqrt(var.counter2))[qn.plot,delta.name]/(cosh(theta.counter2[qn.plot,delta.name]))^2),
            col=c("red","steelblue4","red"), type="l", lty=1, ylim=c(-1,1),
            xlab="Quantile Index", ylab="", main=bquote(rho(y)~"for"~.(delta.para)));
    abline(h=0);
    abline(h=rho.Heckman.counter2, lty=2);
    if(rho.spec==2 & i%%6==0){#mtext(paste("Estimates of Unobserved Selection, ", main.str.counter2, sep=""), outer=T, cex=1.2)
      mtext(main.str.counter2, outer=T, cex=1.2)};
  };
  if(rho.spec!=2){#mtext(paste("Estimates of Unobserved Selection, ", main.str.counter2, sep=""), outer=T, cex=1.2)
    mtext(main.str.counter2, outer=T, cex=1.2)};
}else if(rho.spec==5 | rho.spec==6 | rho.spec==7 | rho.spec==8 | rho.spec==11 | rho.spec==13 | rho.spec==17 | rho.spec==20 | rho.spec==21){
  par(mfrow=mfrow.current,mar=c(5,4,4,2),oma=c(0,0,2,0));
  for(i in 1:length(delta.namelist)){
    delta.para <- delta.namelist[i];   
    matplot(q.plot/100,
            cbind(-lbound.counter2[qn.plot,paste("delta.",delta.para,sep="")],
                  -theta.counter2[qn.plot,paste("delta.",delta.para,sep="")],
                  -ubound.counter2[qn.plot,paste("delta.",delta.para,sep="")]),
            col=c("red","steelblue4","red"), type="l", lty=1, ylim=ylim.current,
            xlab="Quantile Index", ylab="", 
            main=bquote("Effect"~"of"~ .(delta.para)~"on"~delta(y)));
    abline(h=0);
  };
  #mtext(paste("Estimates of Unobserved Selection, ", main.str.counter2, sep=""), outer=T, cex=1.2);
  mtext(main.str.counter2, outer=T, cex=1.2);
}else if(rho.spec==4){
  par(mfrow=mfrow.current,mar=c(5,4,4,2),oma=c(0,0,2,0)); 
  for(i in 1:length(delta.each)){
    delta.para <- delta.each[i];
    delta.name <- paste("delta.",delta.para,sep="");
    matplot(q.plot/100,
            cbind(-tanh(theta.counter2[qn.plot,delta.name])-t(crt.counter2*sqrt(var.counter2))[qn.plot,delta.name]/(cosh(theta.counter2[qn.plot,delta.name]))^2,
                  -tanh(theta.counter2[qn.plot,delta.name]),
                  -tanh(theta.counter2[qn.plot,delta.name])+t(crt.counter2*sqrt(var.counter2))[qn.plot,delta.name]/(cosh(theta.counter2[qn.plot,delta.name]))^2),
            col=c("red","steelblue4","red"), type="l", lty=1, ylim=c(-1,1),
            xlab="Quantile Index", ylab="", main=bquote(rho(x) ~ "for" ~ .(delta.para)));
    abline(h=0);
    abline(h=rho.Heckman.counter2, lty=2);
    if(i%%6==0){#mtext(paste("Estimates of Unobserved Selection, ", main.str.counter2, sep=""), outer=T, cex=1.2)
      mtext(main.str.counter2, outer=T, cex=1.2)};
  };
  if(length(delta.each)%%6 != 0){#mtext(paste("Estimates of Unobserved Selection, ", main.str.counter2, sep=""), outer=T, cex=1.2)
    mtext(main.str.counter2, outer=T, cex=1.2)};
  
  par(mfrow=c(1,1),mar=c(5,4,4,2),oma=c(0,0,0,0)); 
  mycols <- topo.colors(length(rho.name));
  matplot(q.plot/100,
          -tanh(theta.counter2[qn.plot,rho.name]),
          col=mycols, type="l", lty=1, ylim=c(-1,1),
          xlab="Quantile Index", ylab="", main=main.str.counter2);
  #paste("Estimates of Unobserved Selection by Years, ", main.str.counter2, sep="")
  abline(h=0);
  image.plot(legend.only=T, zlim=c(1978,2013), col=mycols, smallplot=c(0.8,0.82,0.2,0.4),
             axis.args=list(at=c(1978,1995,2013),labels=c("1978","1995","2013")));
}else if(rho.spec==15){
  par(mfrow=c(1,1),mar=c(5,4,4,2),oma=c(0,0,0,0));  
  allt <- unique(mysample$year_res);
  mycols <- topo.colors(length(allt));
  matplot(q.plot/100,
          -tanh(replicate(length(allt), theta.counter2[qn.plot,"delta.(Intercept)"]) + 
                  theta.counter2[qn.plot,"delta.year_res"] %o% allt),
          col=mycols, type="l", lty=1, ylim=c(-1,1),
          xlab="Quantile Index", ylab="", main=main.str.counter2);
  #paste("Estimates of Unobserved Selection by Years, ", main.str.counter2, sep="")
  abline(h=0);
  image.plot(legend.only=T, zlim=c(1978,2013), col=mycols, smallplot=c(0.8,0.82,0.2,0.4),
             axis.args=list(at=c(1978,1995,2013),labels=c("1978","1995","2013")));
  
  par(mfrow=mfrow.current,mar=c(5,4,4,2),oma=c(0,0,2,0));
  for(i in 1:length(delta.namelist)){
    delta.para <- delta.namelist[i];   
    matplot(q.plot/100,
            cbind(-lbound.counter2[qn.plot,paste("delta.",delta.para,sep="")],
                  -theta.counter2[qn.plot,paste("delta.",delta.para,sep="")],
                  -ubound.counter2[qn.plot,paste("delta.",delta.para,sep="")]),
            col=c("red","steelblue4","red"), type="l", lty=1, 
            xlab="Quantile Index", ylab="", 
            main=bquote("Effect"~"of"~.(delta.para)~"on"~delta(y)));
    abline(h=0);
  }; 
  #mtext(paste("Estimates of Unobserved Selection, ", main.str.counter2, sep=""), outer=T, cex=1.2);
  mtext(main.str.counter2, outer=T, cex=1.2);
}else if(rho.spec==16){
  par(mfrow=c(1,1),mar=c(5,4,4,2),oma=c(0,0,0,0)); 
  allt <- unique(mysample$year_res);
  mycols <- topo.colors(length(allt));
  matplot(q.plot/100,
          -tanh(replicate(length(allt), theta.counter2[qn.plot,"delta.(Intercept)"]) + 
                  theta.counter2[qn.plot,"delta.year_res"] %o% allt +
                  theta.counter2[qn.plot,"delta.I(year_res^2)"] %o% allt^2),
          col=mycols, type="l", lty=1, ylim=c(-1,1),
          xlab="Quantile Index", ylab="", main=main.str.counter2);
  paste("Estimates of Unobserved Selection by Years, ", main.str.counter2, sep="")
  abline(h=0);
  image.plot(legend.only=T, zlim=c(1978,2013), col=mycols, smallplot=c(0.8,0.82,0.2,0.4),
             axis.args=list(at=c(1978,1995,2013),labels=c("1978","1995","2013")));
  
  par(mfrow=mfrow.current,mar=c(5,4,4,2),oma=c(0,0,2,0));
  for(i in 1:length(delta.namelist)){
    delta.para <- delta.namelist[i];   
    matplot(q.plot/100,
            cbind(-lbound.counter2[qn.plot,paste("delta.",delta.para,sep="")],
                  -theta.counter2[qn.plot,paste("delta.",delta.para,sep="")],
                  -ubound.counter2[qn.plot,paste("delta.",delta.para,sep="")]),
            col=c("red","steelblue4","red"), type="l", lty=1, 
            xlab="Quantile Index", ylab="", 
            main=bquote("Effect"~"of"~.(delta.para)~"on"~delta(y)));
    abline(h=0);
  }; 
  #mtext(paste("Estimates of Unobserved Selection, ", main.str.counter2, sep=""), outer=T, cex=1.2);
  mtext(main.str.counter2, outer=T, cex=1.2)
}else if(rho.spec==18){
  par(mfrow=c(2,1),mar=c(5,4,4,2),oma=c(0,0,2,0)); 
  allt <- unique(mysample$year_res);
  mycols <- topo.colors(length(allt));
  matplot(q.plot/100,
          -tanh(replicate(length(allt), theta.counter2[qn.plot,"delta.(Intercept)"]) + 
                  theta.counter2[qn.plot,"delta.year_res"] %o% allt),
          col=mycols, type="l", lty=1, ylim=c(-1,1),
          xlab="Quantile Index", ylab="", main=bquote(rho(y)~"for Single"));
  abline(h=0);
  #mtext(paste("Estimates of Unobserved Selection by Years, ", main.str.counter2, sep=""), outer=T, cex=1.2);
  mtext(main.str.counter2, outer=T, cex=1.2);
  
  matplot(q.plot/100,
          -tanh(replicate(length(allt), theta.counter2[qn.plot,"delta.(Intercept)"] + theta.counter2[qn.plot,"delta.couple"]) + 
                  (theta.counter2[qn.plot,"delta.year_res"] + theta.counter2[qn.plot,"delta.year_res:couple"]) %o% allt),
          col=mycols, type="l", lty=1, ylim=c(-1,1),
          xlab="Quantile Index", ylab="", main=bquote(rho(y)~"for Married"));
  abline(h=0);
  image.plot(legend.only=T, zlim=c(1978,2013), col=mycols, smallplot=c(0.8,0.82,0.6,0.7),
             axis.args=list(at=c(1978,2013),labels=c("1978","2013")));  
  
  par(mfrow=mfrow.current,mar=c(5,4,4,2),oma=c(0,0,2,0));
  for(i in 1:length(delta.namelist)){
    delta.para <- delta.namelist[i];   
    matplot(q.plot/100,
            cbind(-lbound.counter2[qn.plot,paste("delta.",delta.para,sep="")],
                  -theta.counter2[qn.plot,paste("delta.",delta.para,sep="")],
                  -ubound.counter2[qn.plot,paste("delta.",delta.para,sep="")]),
            col=c("red","steelblue4","red"), type="l", lty=1, 
            xlab="Quantile Index", ylab="", 
            main=bquote("Effect"~"of"~.(delta.para)~"on"~delta(y)));
    abline(h=0);
  }; 
  #mtext(paste("Estimates of Unobserved Selection, ", main.str.counter2, sep=""), outer=T, cex=1.2);
  mtext(main.str.counter2, outer=T, cex=1.2);
};



# Comparison across Counterfactual Levels #
if(rho.spec==1){
  par(mfrow=c(1,1),mar=c(5,4,4,2),oma=c(0,0,2,0));
  matplot(q.plot/100,
          cbind(-tanh(theta.counter1[qn.plot,"delta.(Intercept)"])+t(crt.counter1*sqrt(var.counter1))[qn.plot,"delta.(Intercept)"]/(cosh(theta.counter1[qn.plot,"delta.(Intercept)"]))^2,
                -tanh(theta.counter1[qn.plot,"delta.(Intercept)"]),
                -tanh(theta.counter1[qn.plot,"delta.(Intercept)"])-t(crt.counter1*sqrt(var.counter1))[qn.plot,"delta.(Intercept)"]/(cosh(theta.counter1[qn.plot,"delta.(Intercept)"]))^2,
                -tanh(theta.counter2[qn.plot,"delta.(Intercept)"])+t(crt.counter2*sqrt(var.counter2))[qn.plot,"delta.(Intercept)"]/(cosh(theta.counter2[qn.plot,"delta.(Intercept)"]))^2,
                -tanh(theta.counter2[qn.plot,"delta.(Intercept)"]),
                -tanh(theta.counter2[qn.plot,"delta.(Intercept)"])-t(crt.counter2*sqrt(var.counter2))[qn.plot,"delta.(Intercept)"]/(cosh(theta.counter2[qn.plot,"delta.(Intercept)"]))^2),
          col=c("red","steelblue4","red","red","steelblue4","red"), type="l", lty=c(1,1,1,2,2,2), ylim=c(-1,1), 
          xlab="Quantile Index",ylab="");
  abline(h=0);
  abline(h=rho.Heckman.counter1, lty=2);
  abline(h=rho.Heckman.counter2, lty=2, col="grey");
  mtext(paste("Comparison of Unobserved Selection, ",main.str," (",legend.counter1," vs. ",legend.counter2,")", sep=""), outer=T, cex=1.2);
  
}else if(rho.spec==2 | rho.spec==3 | rho.spec==4 | rho.spec==9 | rho.spec==10 | rho.spec==12 | rho.spec==14){
  par(mfrow=mfrow.current,mar=c(5,4,4,2),oma=c(0,0,2,0)); 
  for(i in 1:length(delta.each)){
    delta.para <- delta.each[i];
    delta.name <- paste("delta.",delta.para,sep="");
    matplot(q.plot/100,
            cbind(-tanh(theta.counter1[qn.plot,delta.name])-t(crt.counter1*sqrt(var.counter1))[qn.plot,delta.name]/(cosh(theta.counter1[qn.plot,delta.name]))^2,
                  -tanh(theta.counter1[qn.plot,delta.name]),
                  -tanh(theta.counter1[qn.plot,delta.name])+t(crt.counter1*sqrt(var.counter1))[qn.plot,delta.name]/(cosh(theta.counter1[qn.plot,delta.name]))^2,
                  -tanh(theta.counter2[qn.plot,delta.name])-t(crt.counter2*sqrt(var.counter2))[qn.plot,delta.name]/(cosh(theta.counter2[qn.plot,delta.name]))^2,
                  -tanh(theta.counter2[qn.plot,delta.name]),
                  -tanh(theta.counter2[qn.plot,delta.name])+t(crt.counter2*sqrt(var.counter2))[qn.plot,delta.name]/(cosh(theta.counter2[qn.plot,delta.name]))^2),
            col=c("red","steelblue4","red","red","steelblue4","red"), type="l", lty=c(1,1,1,2,2,2), ylim=c(-1,1),
            xlab="Quantile Index", ylab="", main=bquote(rho(y) ~ "for" ~ .(delta.para)));
    abline(h=0);
    abline(h=rho.Heckman.counter1, lty=2);
    abline(h=rho.Heckman.counter2, lty=2, col="grey");
    if(rho.spec==2 & i%%6==0){mtext(paste("Comparison of Unobserved Selection, ",main.str," (",legend.counter1," vs. ",legend.counter2,")", sep=""),
                                    outer=T, cex=1.2)};    
  };
  if(rho.spec!=2){mtext(paste("Comparison of Unobserved Selection, ",main.str," (",legend.counter1," vs. ",legend.counter2,")", sep=""),
                        outer=T, cex=1.2)};
}

dev.off()
